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Abstract. In a fractional Cauchy problem, the first time derivative is replaced by 
■ a Caputo fractional derivative of order less than one. If the original Cauchy problem 

governs a Markov process, a non-Markovian time change yields a stochastic solution to 
the fractional Cauchy problem, using the first passage time of a stable subordinator. 
This paper proves that a spectrally negative stable process reflected at its infimum 
has the same one dimensional distributions as the inverse stable subordinator. There- 
£NJ ' fore, this Markov process can also be used as a time change, to produce stochastic 

solutions to fractional Cauchy problems. The proof uses an extension of the D. Andre 
reflection principle. The forward equation of the reflected stable process is established, 
including the appropriate fractional boundary condition, and its transition densities 
are explicitly computed. 



>• . A classical paper of Einstein [13] outlines the link between random walks, Brownian 

motion, and the diffusion equation d t u(x,t) = d^u(x,t). Sokolov and Klafter [33] 
review modern extensions of the theory to encompass random walks with heavy tails, 



O 



1. Introduction 



V) ■ jump processes, and fractional diffusion. This connection between probability and 

differential equations has profound consequences for mathematics [U [7J EJ EH] , and for its 
applications in science and engineering [MJ EHl Elj , including a probabilistic method 
for solving fractional differential equations, by exploiting the associated Markov process 
[37]. More details on fractional calculus, and its connection to probability theory, may 
be found in the recent book of Meerschaert and Sikorskii [26] . 

The time-fractional diffusion equation d^u(x,t) = d^,u(x,t) with < (3 < 1 governs 
the long-time limit of a random walk with power law waiting times ~P(W > t) ~ Ct~^ 
between particle jumps (22J, Theorem 4.2], see also (26j Sec. 4.5]. The limit process 
is a Brownian motion, time-changed by the inverse or hitting time E t of a /3-stable 
subordinator. This time change process is not Markovian, since it remains constant 
over time intervals whose distribution is not exponential. 

The time-fractional diffusion equation is one example of a fractional Cauchy problem. 
More generally, for any strongly continuous Markov semigroup T t f(x) = E x [f(X t )] 
with generator L, it follows from [3j Theorem 3.1] along with (22J Corollary 3.1] that 
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p(x,t) = W[f(XE t )} solves the fractional Cauchy problem 
(1.1) d?p(x,t) = Lp(x,t); p(x,0) = f(x) 



of order < /3 < 1, for any / in the domain of the generator L. The Caputo fractional 
derivative in fll.il) is defined by 



The stochastic process X^ t provides the basis for particle tracking solutions to the 
time-fractional Cauchy problem, in which the underlying process is simulated, and the 
resulting histogram for many sample paths approximates the solution [211 138] - 

The goal of this paper is to provide an alternative method for particle tracking solu- 
tions of fractional Cauchy problems, by establishing an equivalent solution method in 
terms of a Markovian time change. This solves an open problem in [5j Remark 5.2]. We 
will show that a certain Markov process Z t has the same one dimensional distributions 
as the inverse stable subordinator E t , and therefore, one can also compute solutions 
to fractional Cauchy problems by the probabilistic formula p(x,i) = W[f(Xz t )}- The 
proof uses a variation of the D. Andre reflection principle, which may also be useful in 
other contexts. The Markov process Z t is a stable Levy process with index a = 1//3 
and no positive jumps, reflected at the origin. At the end of this paper, we establish the 
forward (Fokker-Planck) equation for this Markov process, which involves a fractional 
boundary condition, and we explicitly compute the transition densities (see Figures [T] 
and . 

Notation. Most of our notation is standard or self-explicatory. Since our paper ad- 
dresses both probabilists and analysts, some words on function spaces and convergence 
seem to be in order. We write CoofO, oo) for the continuous functions which vanish 
at infinity, i.e., lim a; _ il00 u(x) = 0, and this space is always equipped with uniform 
norm ||u|| = sup^o^ |m(x)|. Its topological dual is the space of (signed) Radon mea- 
sures M&[0, oo), and by M ac [0,oo) we mean the absolutely continuous (with respect 
to Lebesgue measure) elements in M&[0,oo). On M^O, oo) we use vague convergence 
which is, topologically, the weak-* convergence; i.e., vaguely if, and only if, 

J udfi n — > fudfi for all u G C^O, oo). The subscripts c,b,ac,oo stand for 'compact 
support', 'bounded', 'absolutely continuous' and 'vanishing at infinity'. 

Fractional integrals and derivatives in the Riemann-Liouville sense are denoted by 



I a and D a , cf. (j4.2|HI4.5p while Caputo derivatives are written as d a , cf. (14. lip . 



In this section, we note a useful extension of the D. Andre reflection principle for 
Brownian motion. Since this extension is not completely standard, we include a simple 



Theorem 2.1 (reflection principle). Suppose that Y t is a Levy process started at the 
origin, with no positive jumps, and let St = sup{K u : ^ u ^ t}. Assume that 
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(1.2) 




2. Reflection principle 



proof. 



P(y > 0) = F(Y 1 ^ 0), for all t > 0. Then 



(2.1) 



P{S t ^x) = P{Y t ^x\Y t ^0) 



P(Y t > x) 
P(Y t > 0) 



for all t, x > 0. 

Proof. Let r x := inf{u > : Y u > x} denote the first passage time process. Since 
(Y t ) t ^o has stationary independent increments, it follows that (Y t+Tx — Y Tx ) t ^ is a Levy 
process, which is independent of the cx-algebra generated by (Yt)t^ Tx , and it has the 
same finite dimensional distributions as (Y t ) t ^ . Consequently 



P (r x ^t,Y t <Y T J = P {r x ^t,Y t - Y Tx < 0) 
= P (t x < t) P (Y 1 < 0) 



(2.2) 




Observe that we have for any cadlag process 



{r x <t}c {S t >x}c {t x ^t}c {S t ^ x} 



for all t and x > 0. Therefore, 



P(S t >x)^ P(r x <: t) 



p(t x < t, y t ^ y Tx ) + p(t x < t, y t < yj 
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On the other hand, 

P(S t > x) ^ F(t x <t) = F(t x < t,Y t ^ Y Tx ) + F(t x < t,Y t < Y Ts 

F(r x <t,Y t ^Y Tx ) 



P (Y 1 ^ 0) 

= p(iTFo) p(T ' <i - y '> K - ) 

= F(Tr>o) P(T - <( - y '> i; " y ''> 0) 

= P(t, <t,Y t >Y Tx \Y, >0) 

>p(r I <t,y l >y T ,|y t >o). 

Therefore, 

P(r x <t,Y t > Y Tx | Y t > 0) < P(5 t > x) < P(r x < t, Y t ^ Y Tx \Y t >0). 
Since Y t has no upward jumps, Y Tx = x and {r x < t} PI {Ft > x} = {y > x}. Therefore, 

p(y > x | y ^ o) ^ p(5 t > x) < p(y ^ x | y ^ o). 

We can now use standard approximation techniques: 

{X ^ x} = f]{X > x - 1/n} = p| {X ^ x - 1/n} 



and 



to get 



{X > x} = |J{X >x+l/n} 

n 

P(y ^ x | y ^ 0) = lim p(y > x - 1/n I Y t ^ 0) 

n— >oo 

^ lim P(5< > x - 1/n) = P(S t ^ x) 



and 



P(S t ^ x) = lim P(S , t > x - 1/n) lim P(y ^ x - 1/n | y > 0) 

n— >oo n— >oo 

= p(y ^ x | y ^ o) 

which proves P(S t > x) = P(y ^ x | y > 0). □ 

Remark 2.2. If (Y t ) t ^ Q is a Brownian motion, then (12. ip becomes the classical reflection 
principle: P(Y t ^ 0) = 1/2, so that Q2H) is equivalent to F(S t > x) = 2P(y ^ x). 

The proof of Theorem 12 . 1 1 relies essentially on local symmetry and the strong Markov 
property. Let y be a strong Markov process with cadlag paths and transition function 
Pt(z, dy) = P 2 (y G dy), and write r 2 = inf{w > : Y u — z > x} for the first passage 
time above the level x + z for the process Y t started at z; observe that, in general, 
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Yt* x + z. We can use the strong Markov property in ( 12. 2ft to get for any starting 
point z 

P 2 {r z ^ t,Y t < Y T *) = P z (r* ^ t,Y t -Y T * < 0) 



/ p%M (y (w) - y < 0) P 2 (<L;). 
./fr*<ti 



'{*2<t} 

If we assume, in addition, some local 'symmetry', i.e., that for some constant c G (0, oo) 
we have 

(2.3) F2 j Ft - - < °j = c for all t > 0, z G R, 

V ; P z (Y t -z^O) 

then we get P Y -lH (Y t _ T * {w) - Y < 0) = c P y ^ } (y_ r , H - F ^ 0) and, with a similar 
argument, 

P 2 (t z x ^t,Y t < Y TS ) = cP 2 (r| ^ t, y ^ y.) . 

This means that we can follow the lines of the proof of Theorem 12.11 to derive the 
following general result. 

Theorem 2.3 (Markov reflection principle). Suppose (y,P z ) is a strong Markov pro- 
cess satisfying the local symmetry condition (12. 3p . Set S t = sup{y u — Y : ^ u ^ t}. 
Then we have for all t, x > and zGR 

p 2 (^ > x) ^ P z (s t > x, y ^ y r| | y ^ z ) 

«S P 2 (^i -z ^x\Y t ^z). 

IfY t has only non-positive jumps, then Y T z = x + z a.s., and we get for all t, x > and 
zeR 

(2.5) F z (S t ^x) =P z (Y t - z^ x\Y t > z). 

3. An equivalent Markov time change 

Suppose that D t is a standard stable subordinator, a Levy process with E[e _sDi ] = 
e -tsP f Qr a jj t ^ 0, and define its inverse (hitting time, first passage time) process by 

(3.1) E t = inf{r > : D r > t} 

for all t ^ 0. In this section, we prove that a certain Markov process Z t has the 
same one dimensional distributions as the inverse stable subordinator (13. ip with index 
(3 G [1/2, 1). Then we use the Markov process Z t as a time change, to develop efficient 
probabilistic solutions to fractional Cauchy problems. 

To construct this Markov process, first consider a stable Levy process Y t with char- 
acteristic function 

(3.2) E[e ifcy <] = e t{lk)a 

where a = 1/(3. If (3 = 1/2, then a = 2, and Y t is a Brownian motion with variance It. 
Now define 

(3.3) Z t = Y t - inf{y : ^ s < t}. 
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The reflected stable process (13.31) is also the recurrent extension of the process Y t killed 
at zero, which instantaneously and continuously leaves zero, see Patie and Simon [29] . 
Let Coo(R) denote the Banach space of continuous functions / : R — > R that tend to 
zero as \x\ — > oo, with the supremum norm. We say that a time-homogeneous Markov 
process X t is a Feller process if the semigroup T t f(x) = F l [f(X t+s )\X s = x] satisfies 
T t f £ Coo(R) and T t f — > f as t — > in the Banach space (supremum) norm, for all 

/ e CUR). 

Lemma 3.1. Let E t be the inverse stable subordinator (13. ip . where D t is a standard 
(3 -stable subordinator with E[e~ sDt ] = e~ tsf> and (3 £ [1/2, 1). Define Z t by ( I3.3p . where 
Y t is a stable Levy process with characteristic function (13. 2p . with index a = 1/(3. Then 
Z t is a Feller process, with the same (one- dimensional) marginal distributions (given 
Zq = 0) as E t , for any fixed t ^ 0. 

Proof. Define the running infimum I t = inf{Y, : ^ s ^ t} and the running supremum 
St = sup{F s : ^ s ^ t}. Let Y t = —Y t denote the dual process, and let I t and 
St denote the running infimum and supremum of Y t , respectively. Since Y t is also a 
Levy process, it follows from [TU| Section VI. 1, Proposition 1] that St — Y t is a Feller 
process, and since S t = —I t , it follows that S t — Y t — —It + Y t — Z t . Hence Z t is a 
Feller process. Since the standard stable subordinator D t has characteristic function 
gj-gifcDtj _ Qxp^—tl—ik) 1 ), an application of the Zolotarev duality for stable densities 
Theorem 4.1] implies that 

(3.4) P(E t > x) = F(Y t > x\Y t ^ 0) for all t > and all x > 0. 

The stable Levy process Y t has no positive jumps, and since Y t is self-similar (e.g., see 
[261 P- 105]), it follows that F(Y t > 0) = Vit 1 '^ > 0) = ¥{Y 1 > 0) for all t > 0. 
Since every stable process has a Lebesgue density (e.g., see (261 P- 107]), we also have 
P(Yf = 0) = for all t > 0. Then it follows from Theorem 12 . 1 1 1 hat 

(3.5) F(Y t > x\Y t ^ 0) = F(S t ^ x) for aU t > and all x > 0. 
Finally, a fluctuation identity [TUl Section VI. 1, Prop. 3] implies that 

(3.6) F(S t ^ x) = F(Z t > x) for all t > and all x > 0. 

Then the theorem follows by combining (13.4)) . (13. 5p . and (I3.6p . □ 

Remark 3.2. It follows from the Zolotarev duality formula for stable laws [201 [39] l na l 
in fact F(Y t ^ 0) = 1/a, see Theorem 4.1]. 

The next result shows that fractional Cauchy problems can be solved by a proba- 
bilistic method, using a Markovian time change. 

Theorem 3.3. For any f3 £ [1/2,1), let Z t be given by ( 13 .3|) . where Y t is an a-stable 
Levy process with characteristic function (13.2)) for a = 1/ (3. If X t is an independent 
Markov process such that T t f(x) = ~E x [f(X t )] forms a uniformly bounded, strongly 
continuous semigroup on some Banach space B of real valued functions, with generator 
L, then p(x,t) = Ei x [f(Xz t )} solves the fractional Cauchy problem ( II. ip for any f £ 
D(L), the domain of the generator. 
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Proof. A general result [31 Theorem 3.1] shows that, if L is the generator of a uniformly 
bounded, strongly continuous semigroup on some Banach space B, and if u(-,t) G B 
solves the Cauchy problem d t u(x,t) = Lu(x,t); u(x, 0) = f(x) on that space for some 
/ G Dom(L), then the solution to the associated fractional Cauchy problem (11. ip on 
that Banach space is given by 



(3.7) p(x, t) — / u(x,r)h(r,t) dr 

Jo 

where 

(3.8) h(r,t) = tr- 1 - 1 "* g p (tr-W) 
and gp(t) is the function with Laplace transform 

POO 

£g(s)= / e~ st g p {t)dt = e- sP 
Jo 

for some < < 1. Of course gp(t) is the probability density function of a standard 
/3-stable subordinator. Let D t be the associated Levy process with E[e _sDt ] = e~ tsP 
and apply [221 Corollary 3.1] to see that (I3.8P is also the probability density of the hit- 
ting time (13. ip . Since T t f(x) = W[f(X t )} is a uniformly bounded, strongly continuous 
semigroup, it follows that p(x,t) = E x [f(XE t )} solves the fractional Cauchy problem 
(11. ip . Lemma [3.11 implies that the random variables E t and Z t have the same distribu- 
tion for any t > 0, and we also have Z = E = almost surely. Then it follows that 
W[f{X Zt )} = E x [f(X Et )}, which completes the proof. □ 

Remark 3.4. Theorem 13.31 confirms a conjecture in the paper [3 Remark 5.2]. There 
we set Z t = Fcr(t) where a(t) = mi{u > : H u > t} and H u = J ly a>0 c?s. In essence, 
the negative excursions are cut away, and the positive excursions are joined together 
without any gaps in time. Since Y t has no positive jumps, any up-crossing at the origin 
is a renewal point, so this process has the same distribution as ( 13. 3p . 



Remark 3.5. Since the stable Levy process Y t with index a G (1,2] and characteristic 
function ( 13. 2 p has no positive jumps, the first passage time D t := inf{r > 0; Y r > t} 
is a stable subordinator with E[e~ sDt ] = exp(— ts 1//Q ), and the supremum process E t = 
sup{F r : ^ r ^ t] is the inverse stable process (13. ip of D t , see Bingham [12]. This 
allows us to compare the processes Z t and E t in terms of their sample paths. It can 
also be used to give an alternative proof of Lemma 13.11 Since S t = E t , we certainly 
have P(-Et ^ x) — P(5j ^ x), and then the result follows by the fluctuation identity 
(13. 6p . For the case a = 2, Z t is the reflected Brownian motion, and E t is the supremum 
of that same Brownian motion. 



4. The reflected stable process 

The equivalence of one dimensional distributions established in Lemma 13.11 also has 
some interesting implications regarding the reflected stable process Z t defined in ( 13. 3p . 
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Proposition 4.1. Let Y t be a stable Levy process with characteristic function ( 13. 21) 
and index a G (1,2). Then the reflected process Z t in (13.31) has the right- continuous 
probability density 

( t^x-^^g^tx- 1 ^) x > 

(4.i) P (x, t) = I r^/r(i -p) x = o 

[0 x < 

for any t > 0, where /3 — 1/a and gp(f) is the probability density function with Laplace 
transform Lg(s) = exp(— s 1 /"). 

Proof. Lemma 13.11 shows that the probability distribution of Z t is also the probability 
distribution of the random variable E t , the hitting time (13.11) of a stable subordinator 
D t such that E[e" sDt ] = exp(-ts^), and = 1/a. Then [221 Corollary 3.1] implies 
that E t has a Lebesgue density h(x,t) = tfi~ x x~ x ~ x 't ] g^tx -1 '^) for x > 0. Since D t 
is right-continuous, it follows from (13. ip that E t > for t > 0, so that p(x,t) = for 
x < and t > 0. Finally, since g^(x) ~ /3x _/3 ~ 1 /r(l — f3) as a: — )■ oo, it is easy to check 
that p(x, t) ->• t~ p /T(l - (3) as x ->■ 0+ for any i > 0. □ 

The next result shows that the density (14.11) of Z t started at Z = solves a fractional 
boundary value problem. Given a real number a > which is not an integer, define 
the positive Riemann-Liouville fractional integral 



(4.2) i:f\x) = — / /(x - y)^" 1 dy=— f(y)(x - y)^ 1 dy, 

1 (a) 7o 1 («) J-oo 

the negative Riemann-Liouville fractional integral 

(4.3) I«J( X )= f{x + y) y^ dy= f{y){y-xY- l dy, 

Y[a) J T(a) J x 

the positive Riemann-Liouville fractional derivative 

d n 1 d n r x 

(4.4) D^m:=—ir a f(x) = - r— / mtx-yr-^dy, 

dx n 1 (n — a) dx n J _ 00 

and the negative Riemann-Liouville fractional derivative 

< 45 > == - f& £ r - *• 

where n — 1 < a < n. If a G (1, 2), then n = 2. 

Theorem 4.2. Lei It fre a stable Levy process with characteristic function (13 .2D and 
index 1 < a < 2, and /ei Z 4 &e ine reflected process CI3 .31) . Tnen i/ie probability density 
(14. ip o/ Zt solves the fractional boundary value problem: 



(4-6) 



d t p(x, t) = D a _ x p{x ) t) for t > and x > 0; 

= Dl-t-pfa t) fort>0 and x = 0. 
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Proof. Let q(x,t) denote the probability density of the stable Levy process Y t . An 
application of the Zolotarev duality formula for stable densities [5] Theorem 4.1] to- 
gether with Lemma 13.11 shows that (3p(x,t) = q(x,t) for all x > and t > 0, where 
f} = F(Y t > 0) for any t > 0. It follows from [261 Example 3.29] that the stable density 
q (x, t) solves 

(4.7) d t q{x,t) = D« x q{x,t) 

for all t > and all x G R. Hence we also have d t p(x, t) = D" x p(x, t) for t > and 
x > 0. When x = 0, it follows from Proposition 14.11 and the continuity of x H- q{x,t) 
that /3p(0,t) = q(0,t) for t > 0. Then d t p(x,t) = D" x p(x,t) at x = as well, since 
D" x p(x, t) at x = depends only on the values of p(x, t) for x ^ 0. 

It remains to verify the fractional boundary condition. From the general definition 
(14 .5p it follows that for 1 < a < 2 we have 

1 d 2 r°° 

(48) D "*^w^)^L ,W(! - ir * 

Since 1 = ~P{Z t ^ 0) is constant for all t > 0, it follows that 

POO 

= d t F(Z t ^ 0) = / d t p(x, t) dx 
Jo 

D°: x p(x, t) dx 

a r°° \ 
P(y,t)(y - x)^ a dy J 

= -y -^-DTM^t) dx 

since lim^oo D°L~ p(x, £) = alim^oo D"~ 1 q(x, t) = 0, by an elementary estimate using 
the fact that q(x, t) ~ Cx~ a ~ l for some C > as x — > oo. □ 

Let be the stochastic process defined in (13.31) . where Y t is a stable Levy process 
with index 1 < a < 2 and characteristic function (13. 2p . Since Z t ^ by definition, it 
follows from Lemma 13.11 that Z t is a conservative time-homogeneous Markov process 
whose (backward) semigroup T t f(x) = Wi[f(Z t+s )\Z s = x] is strongly continuous and 
contractive (i.e., \\T t f\\ ^ ||/||) on the Banach space X = C^O, oo). Let A denote the 
generator of this semigroup, with domain D(A). Next we will show that this semigroup 
is also analytic and give a core for the generator. Recall that a core Ca of a closed linear 
operator A is a subset of its domain that is dense within the domain if the domain is 
equipped with the graph norm; i.e., Ca is a core of A if for each / e D(A) there exists 
a sequence {f n } C CU such that / n — )■ / and A/ n — >■ A/. 
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(4.9) _ P d / 1 d 

dx \r(2 — a) dx 



Write 

S b := {/ £ ^[0, oo) : /" e C(0, oo), /"(x) = 0(1) as x -> oo, 

/" (x) = 0(x a ~ 2 ) as x -)• 0, /' £ C 6 (0, oo), /' (0+) = o} 



(4.10) 



and denote by the (positive) Caputo fractional derivative of order a > 0, which can 
be defined by 

(4.11) d a J{x) = ir a f {n \x) = 1 f\x - y) n - l - a f {n \y) dy 

T(n - a) J 

where n— 1 < a < n and is the nth derivative of /. The Caputo fractional derivative 
differs from the Riemann-Liouville form f)4.4p because the operations of differentiation 
and (fractional) integration do not commute in general. For example, when < a < 1 
we have d%f(x) = D%f(x) — /(0)x~ a /T(l — a) for suitably nice bounded functions (e.g., 
see [261 P- 39]). 

Theorem 4.3. Let Z t denote the Feller process (13. 3p where Y t is a stable Levy pro- 
cess with index a = 1/(3 £ (1,2) and characteristic function (13.21) . and let T t f(x) = 
Ei[f(Z t+s )\Z s = x] be the associated transition semigroup on C oo [0,oo) with generator 
A. Then {T t } t ^ Q is analytic, Ca = {/ £ S& : d%f £ 0^0, oo)} is a core of A and 
Af = d£f for all f £ Ca, where is the Caputo fractional derivative (14.111) . 

Proof. It follows from Proposition 4] that C A C D(A) and Af = d%f for all / £ C A - 
Note that the extension from /" bounded to |/"(x)| = 0(x a ~ 2 ) at x — 0+ is also 
mentioned in the proof. 

Since A generates a strongly continuous contraction semigroup, the resolvent opera- 
tors R(X, A) := (A — A)~ l exist for all A with Re A > and they are bounded operators. 

Let S e be the span of {x (-> exp(— cx) : c > 0}, which is dense in CoofO, oo) by a 
standard Stone- Weierstrafi argument. If we can show that R(X,A)S e C Ca for some 
A > 0, then Ca is a core of A. 

Indeed, pick / £ D(A) and g = \ f — Af. Since S e is dense in Coo[0, oo), there exists 
a sequence {g n } C S e with g n ->■ g. Thus, /„ = R(X, A)g n -> / and Af n = \ f n - g n -> 
\f — g = Af. Since / n £ i?(A, A)^ e C Ca, we see that Ca is a core. 

It remains to show that f\ = R(X, A)g is an element of Ca for any g £ S e . The Laplace 
transform of a bounded and measurable function g is the unique analytic function 
defined by Lg(s) = J °° e~ sx g(x) dx for Res > 0; in particular, if g(x) = e~ cx for some 
c > 0, its Laplace transform is given by £g(s) = l/(s + c) for any s with Re s > 0. 

We will first define the function fx via its Laplace transform, and later we will show 
that in fact fx = R(X, A)g. Let 

1 

(4.12) f x {8) := ^ " ^ + C 



A - s a 

Next we show that fx is the Laplace transform of a function /a £ Sb- Note that /a can 
be extended to an analytic function on C\ (— oo,0], since the singularity at s = A 1 /" 
is removable. Furthermore, sf\(s) is bounded on the complement of any neighborhood 
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of s = — c, and hence fx is the Laplace transform of a bounded analytic function f\ 
on (0,oo) in view of [2j Theorem 2.6.1] with u = 0. Using the Tauberian theorem 
[21 Theorem 2.6.4] we obtain that hm^^ f\(x) = linx^o sLf\(s) = 0, \im x ^ f\(x) = 
Hindoo sLf\(s) = A 1 / Q_1 /(A 1 / Q + c), and hence fx G C^O, 00). Furthermore 

x \ 1/a— 1 , e 



(4.13) £ f' {s) =sL /a(s) _ /a(Q) = s + c A^+c _ = s + c Ai/°+c 

Using the Tauberian theorem again, / A (0) = and f' x G Coo[0, 00). Taking the Laplace 
transform of f 2 := f' x ' + A i^ +cr (ali) y ields 

« Jif s\*-/ a nc l— a na 2 na \l/et r-a^~ a 

p hi/ \ . Cb _ s+c A'/°+c , Cb _ C6 — CSA Ch 



AV" + C A-s Q AV« + c (s + c^A 1 /" + c)(A — s a ) \ l l a + c 

and then the Tauberian theorem implies that linx^oo / 2 (aO = 0, lim^o^^) = 0. 
Hence fx G S b . 

Next we show that fx G Ca- Since the function x H- x p l[o !00 )(a ; ) has Laplace trans- 
form s~ p ~ l Y(p + 1) for any p > —1, it follows from ( 14. lip along with the convolution 
theorem for the Laplace transform that for / G Sj, we have 

= (s 2 £/(s) - sf(0) - f (0)) = S ^/( S ) - S a "V(0). 

Taking the Laplace transform of \f\ — d%f\, we therefore obtain 

XLfx(s) - £[d:fx](s) = XLfx(s) - s a Lf x (s) + s^f^O) 

1 s a-l X l/a-l ^ i A l/ a -l 

( 4 - 14 ) ~7T^~ AV- + C +S AV« + C 

1 

i 5 
S + C 

which is the Laplace transform of g. Since /a and g are elements of C^O, 00), we 
see that <9"/a G Coo[0, 00) and hence fx G Ca- Now it follows from the uniqueness of 
Laplace transforms that fx = R(X,A)g. Therefore R(X,A)S e C Ca and hence Ca is a 
core. 

Finally we show that {T t }t^o is an analytic semigroup. Recall that the resolvent 
R(X, A) is a bounded operator for all Re A > 0. Then a general result from the theory 
of semigroups (21 Corollary 3.7.12] states that {T t } is bounded (i.e., for some M > 
we have ||T(£)|| ^ M in the operator norm for all t ^ 0) and analytic if ||Ai2(A, A)\\ is 
bounded in the right half plane; i.e. if there exists some M > such that 

(4.15) \\XR(X,A)g\\^M\\g\\ 

for all Re A > and all g G Coo[0, 00), recalling that \\g\\ = sup^Q \g(x)\ is the Banach 
space norm. Note that it is actually enough to show ( 14. 15}) for g in a dense subset (such 
as S e ) due to the continuity of R(X, A). 

Let g G S e , and recall that the Laplace transform 



POO 

(4.16) L[R(\,A)g](s)= e~ sx R(X, A)g(x) dx 

Jo 



A 
11 



can be extended analytically (and uniquely) to s G C \ (—00, 0] for any Re A > 0. Then 
for x > 0, the complex inversion formula for Laplace transforms [Ml Theorem II. 7.4] 
implies that 

1 f a+iN Lg(s) - if^LgiX 1 / ) 

(4.17) [R{\,A)g](x) = —tim / e" 9K ' ^ 9K Us 

2m N^ooJ a _ iN A-s a 

independent of a > 0. For g G S e , the integrand in (14.171) is bounded and analytic on 
an open neighbourhood of the set {s G C : ^ Re(s) ^ a, \s\ > r} for any < r < a. 
Then a standard argument using the Cauchy integral theorem shows that (I4.17P also 
holds for a = 0. Substituting z£ = s, we have 



2vr n-k*> J^ N I A - (z£) a A - (i£) c 
(4.18) = J_ lim f e ^ f MO _ ftp"- 1 ^(A^) \ 

= — lim / e «*ji(£) # _ J_ lim / e^J?(g) dC ^rrTT^ 

2-KN^oo _ N XK ' 2TlN^oo _ M XK ^ J ^ A 1 " 1 /" 



for any Re A > 0. Observe that 



is the Fourier-Laplace transform of the marginal densities of an a-stable process, where 
J e~ l ^ x g a (x) dx = exp((i£) Q ). They are also the Green's functions (convolution kernel) 
for the bounded analytic convolution semigroup generated by D" on Li(R) and also, 
by the transference principle in (U Theorem 4.6 and Corollary 4.2], on Coo(R). Hence 
the resolvent (the Laplace transform of the semigroup) satisfies 



1 f°° 

R{\D«)g{x) = - e^I l x (0d£ 

2n J -OO 



and then [23 Corollary 3.7.12] implies that \\\R(\,D%)g\\ ^ M\\g\\ for some M > 0, for 
all Re(A) > 0. 

In order to estimate the term involving If, observe that, since g a is an a-stable 
density, g a (x) is differentiable, g a {x) ~ x~ a ~ l as x — > 00 and g a (x) decays super- 
exponentially as x — > — 00; hence the function 



xt x 1 fx 



is Lebesgue-integrable in x and differentiable for x > 0. Use (I4.19p to see that 
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Since /f (£) is the Fourier transform of an integrable function which is also different iable, 
the Fourier inversion theorem [361 Theorem7.2-l] yields 



— lim 

2n Af-i>oo 



A 1 



■F x (x) 



N 



a 



for all x > 0. Since |if(£)| ~ \C\ 1 as |£| ~~ * °°; one needs the Cauchy principal value 
here. A simple change of variable then yields 



—F x (x) 
a 



-M 



X 



X 



dt 



e~ w g a (u) du 



C 1 



for all Re (A) > and x > 0. Since |A 1/o | < Re(A 1 /")/ cos(vr/2a) for Re(A) > 0, and 
Re(A|£/i(A)|) ^ Halloo for any h G Coo[0,oo), it follows that 



X- 



A l-l/a 

and hence 



lA^^A 1 



cos(7r/2o;) 



Re(A 1/a )£#(A 



cos(7r/2a) 



||Ai2(A,A)p|| < M 



cos(7r/2a;) 

for all Re A > 0. Then it follows from [2] Corollary 3.7.12] that {T t }t^o is analytic. □ 

Remark 4.4. Patie and Simon [29] show that the reflected stable process Z t in Theo- 
rem 14.31 has the backward generator 



(4.20) 



Af(x) = /'(0) 



l-Q 



rf2-a 



+ 



yl-a 

f'(x-y)— -dy. 

o r(2-a) 



They also give the exact domain of the generator [29l Proposition 2.2]. If / G Sb, then 
/'(0) = 0, and Af reduces to the Caputo fractional derivative (14.11]) . 



In view of Theorem 14 .3[ T t f(x) = ~E[f(Z t+s ) \Z S = x] is a strongly continuous, analytic 
semigroup on the Banach space X := Coo[0, oo) with the supremum norm, with gener- 
ator Af(x) = d%f(x) for / G Sf, such that d%f G X. The dual (or adjoint) semigroup 
T t * is defined on the dual space C*[0, oo) = M&[0, oo) of finite signed Radon measures 
on [0, oo) equipped with the total variation norm: Given a measure fi G M&[0, oo), use 
the Jordan decomposition to write \x = fi + — pT uniquely as a difference of two positive 
measures, and define = /i + [0, oo) + /i~[0, oo). The dual semigroup satisfies 



(4.21) 



T t f(x)n(dx) 



f{x)[T^}{dx) 



for all / G Coo[0,oo) and all \i G M&[0, oo). See [HI Section 2.5] for more details. 
In probabilistic terms, since T t f(x) = ~E[f(Z t + s )\Z s = x] for this time- homogeneous 
Markov process, equation (I4.2ip implies that 



T t f(x)fi{dx) 



E[f(Z t+s )\Z s = x]fi(dx) 
f(y)P t (dy,fi)= [ f(y)[T t *fi](dy) 
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where Pt(y,ji) = J P(y,x,t)fi(dx) and P(y,x,t) = F[Z t+s ^ y\Z s = x] is the transition 
probability distribution of the Markov process Z t . Hence, if /i is the probability distri- 
bution of Z s , then T^fi(dy) = Pt(dy, /i) is the probability distribution of Z t+S . The dual 
semigroup is also called the forward semigroup associated with the Markov process Z, 
since it maps the probability distribution forward in time. 

Next we will compute the generator A* of the forward semigroup. This is the adjoint 
of the generator A of the backward semigroup, in the sense that 

J Af(x)fx(dx) = J f{x)[A*n\{dx) 

for all / G D(A) and fi G D(A*). Theorem 14.51 will show that every measure fi G D(A*) 
has a Lebesgue density g G L^O, oo), so that fi(dy) = g(y)dy, and that the adjoint 
A*g := A* n of the positive fractional Caputo derivative Af(x) = d"f(x) in our setting 
is the negative Riemann-Liouville fractional derivative A*g(y) = D°_ y g(y) using (I4.5p . 

The forward semigroup T t * of a Markov process is not, in general, strongly continuous 
on Mb[0, oo). That is, there exist measures /i such that T t */i -/> /i in the total variation 
norm as t J, 0. For example, if T t * is the forward semigroup associated with the diffusion 
equation d t p = d^p, and /i = 5 is a point mass at the origin, then T 4 */i is a Gaussian 
probability measure with mean and variance 2t for all t > 0, and since yu{0} = 1 and 
T 4 /i{0} = for all t > 0, we have \\T t fi — =1 for all t > in the total variation 
norm. 

To handle this situation, we introduce the sun dual space of X := Coo[0, oo): 

X Q := {/i G X* : lim \\T*a - /ill = 0} 

is closed subspace of X* = M&[0, oo) on which the forward semigroup is strongly con- 
tinuous. It follows from basic semigroup theory [TU Section 2.6] that for fi G X , 
T t *fi G X® for all t ^ 0, and X = D(A*). The restriction of {T*} t>0 to X is called 
the sun dual semigroup {T®}t^o with generator A & fi = A*n for all \i G D(A & ), where 

(4.22) D(A Q ) ={/i£ D(A*) : A*fi G X }. 

For the reflected stable process, we will show in Theorem 14. 51 that C [0, oo) is the space 
of absolutely continuous elements of M&[0, oo), 

M ac [0,oo) = {n G M 6 [0,oo) : fi(dy) = g(y) dy, 3g G ^[0, oo)}. 

and we will derive the forward equation of the reflected stable process on the sun-dual 
space. For a general bounded measure \i G M&[0, oo), we will then prove in Corollary 14. 81 
that T*fi can be computed as the vague limit of T /i n , where fi n — > ji vaguely, and 
fi n G C [0, oo) for all n. 

Theorem 4.5. Let Z t denote the Feller process ( 13 .3|) . where Y t is a stable Levy process 
with index a = 1/(3 G (1,2) and characteristic function ( 13. 21) . with (backward) semi- 
group T t f(x) = Fj[f(Z t+s )\Z s = x] on Coo[0, oo). Then C [0, oo) = M ac [0, oo) and the 
generator A Q g := A & fi of the sun-dual semigroup {T°(t)} t>0 is given by 

(4.23) A Q g(y) = D a _ y g{y) 

with domain D(A Q ) = {g G L^oo) : D* y g(y) G L^O, oo), D^g^) = 0}. 
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Proof. Suppose A*ji = v G M&[0,oo) for some yU G Mb[0,oo), so that / Af{x)\x{dx) = 
J f(x)u(dx) for all / G D(A). Set v(x) := i/[0,x] for x ^ and v(x) = for x < 0. 
If / G S b with 9£/ G CoofO, oo), then it follows from Theorem that / G D(A) and 
Af = d%f. It is obvious from the Definition (T4TTT]) that d™' 1 f (x) = d£f(x). Let 

Jo 1 (") 

denote the positive Riemann-Liouville fractional integral ( 14. 2 p of order a > for a 
function / G Coo[0, oo), and apply the general formula [El Eq. (1.21)] 7" _1 9" _1 / / (x) = 
f'(x) - f (0) to see that f(x) = I^d"' 1 f'(x) = I%- l d«f{x) = I^ x Af{x). Since 
f(x) —> as x — > oo, and u(0) = for x < 0, we can apply the integration by parts 
formula fT7\ Theorem 19.3.13] 

1 ' f{x)v{dx) = f{b)v{b) - f{a)v{a) - f v{x)f{x)dx 



with a < 0, and then let b — > oo, to see that 

POO 

f{x)v(dx) = — I f'(x)v(x)dx. 
o Jo 

Thus, for all / G S& with <9"/ G Coo[0, oo), a Fubini argument yields 

/.oo /-X / _ \a-2 

/(x)i/(dr) = - / / 1 ^ A/(y) dyt^) dx 
Jo Jo r (« - 1) 

/•oo /-oo / _ \a-2 

(4.24) = - / / 1 y; v(s) dx A/(y) rfy 

Jo Jy T ( a - !) 



oo 



o 



Next we will show that S 1 := {Af : f G Sj with d£f G Coo[0, oo)} is dense in 
Coo[0, oo), and then it will follow that any measure /x G D(A*) has a Lebesgue density 

(4.25) <?(</) = -J l r( tt - l) <g) ^ = 

where = z/ and t> (x) = i/[0, a;]. Let C£°[0, oo) denote the space of smooth functions 
with compact support, i.e., such that h(x) = for all x > M, for some M > 0. It is 
not hard to check that the space Q = {h G C^°[0, oo) : J /i = 0} is dense in Coo[0, oo). 
Then we certainly have 

r M ( r _ ^a-l _ T a-1 

lim I£h{x) = lim / ^ '—— h(s) ds = 0, 

x->oo x->oo J T(a) 

and therefore, the function f(x) = I™h(x) is an element of CoofO, oo) for any h G Q. 
Elementary estimates suffice to check that / G S& as well. Since the positive Caputo 
derivative is a left inverse of the positive Riemann-Liouville integral (6J Eq. (1.21)], we 
also have d£f(x) = d*I«h(x) = h(x) G C^oo). Hence h = Af G S for all h e Q, 
and thus Q C S. Now for any / G CoofO, oo) there exists a sequence / n — > f in the 
supremum norm, with f n ES for all n. Then some simple estimates can be used to 
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verify that J f(y)fi(dy) = J f(y)g(y)dy, and it follows that the measure \i has the 
Lebesgue density g. Hence we can identify D(A*) with a subspace of L^O, oo). 
Next we will show that this subspace is dense in L 1 [0, oo). Define 

■= rfey (J - ^T' 1 l[o,i/n)(a:), 

and note that (f) n (x) = 1 for x E (0, 1/n) by a straightforward computation. Note 
that the set U := {g(x) - [D^- 1 g{Q)]^ n : g E G c °°[0, oo), n E N} is dense in L^O.oo), 
and that D^ x g(x) G L x [0, oo) for all g G U. Furthermore, U is a subset of 

G={geL x [0,oo) : D" x g G L 1 ^, oo), D»;^(0) = 0}. 

For any g E G and f E Sb with <9£/ G G^O, oo), we have D^^(0) = and f'(0) = 0. 
Then Theorem I4.3[ a Fubini argument, and integration by parts (twice) using equation 
ffl~5j) yields 

y {y-xf- a 



Af(y)g(y)dy= / / v * f"(x) dx g(y) dy 

to Jo r(2 - a) 

' a g(y)dyf»(x)dx 











r(2-a 

(4-26) =/ I^g(x)f(x)dx 

DT x 1 g(x)f'(x)dx 
f(x) D°_ x g{x)dx. 



o 



oo 







Furthermore, as {/ 6 Sj, : A/ G Goo[0, oo)} is a core, for all / G D(A) there exists a 
sequence f n in the core such that f n —>f and A/ n — > Af. Hence (I4.26P holds for all 
/ G D(A) and therefore (7 G D(A*) and A*g = D°_ x g for any g E G. Since C/ is dense 
in ^[0,00), and U C G C D(A*), it follows that D(A*) is dense in ^[0,00). Since 
G®[0, 00) is the smallest closed set containing D(A*) by definition, and since L 1 [0,oo) 
is a closed subspace of M&[0, 00), we have shown that L 1 ^, 00) = G® [0, 00). 

Equation (14.221) implies that ^ = A*fi is an element of G® [0, 00) for any /x G D(yl ), 
and therefore, we have z/(dx) = h[x)dx for some /i G L 1 [0,oo), as well as fi(dy) = 
g(y) dy for some g E L 1 ^, 00). Since v(x) = u[0,x], it follows that h{x) = v'{x). Since 
D®' 1 is a left inverse of J"" 1 in general, it follows from (I4.25P that v(y) = —D a ^ y 1 g{y). 
Then we have 

(4.27) h(x) = v'(x) = A [-£>^(x)] = ^(x) = A*^(x) 

for all g E D(A Q ), which proves the generator formula f!4.23j) . Equation (14.271) also 
shows that D" x g(x) E L^C^oo) for all g E D(A G ), and since v is continuous with 
v(0) = 0, it follows that D2~^(0) = v(Q) = for all g E D(A Q ). This proves that 
D(A Q ) C G. Since D(A ) is defined as the set of g E D(A*) such that A*g E L l [0, 00), 
it follows from fl426|) that g E D(A & ) for all g E G, so that G C D(A ) as well, which 
completes the proof. □ 
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Theorem 14.51 establishes the forward equation of the reflected stable process Zf, for 
certain initial conditions. It shows that, for any initial condition Ho{dx) = g(x)dx where 
g G L l [Q, oo), /it '■= T Q fi solves the Cauchy problem 

d tf i(t) =A /i(t); t M(0)=fi 

on the sun-dual space. This implies that, for any probability density po(x) such that 
Po(%) — for x < 0, the function p(x,t) = T e po(x) solves the forward equation 

(4.28) d t p(x,t) = D« x p(x,t); p(x, 0) = p (x), D°^p{x,t) = 0. 

x=0 

Remark 4.6. Proposition 14.11 implies that the probability distribution of Z t started at 
Z = has a density function p(x, t) given by (14. ip . Fix s > and take po(x) = p(x, s). 
Then (14.281) holds for all t > s. Since s > is arbitrary, this can be used to give an 
alternative proof of Theorem 14.21 

The next two results will allow us to compute the transition probability density 
y i — y p(x, y, t) of the time-homogeneous Markov process y = Z t+S for any initial state 
x = Z s , by applying Theorem 14.51 to a sequence of initial conditions fi n G X G such that 
A*n ~~ >* $x- The first result shows that the transition density exists. 

Corollary 4.7. The semigroup T® is a strongly continuous bounded analytic semigroup 
on L 1 (R + ) and the transition probability distributions T£5 X have smooth densities y H- 
p(x, y, t) for all t > and all x ^ 0. 

Proof. It is well-known that the spectra of A and A* coincide, and R(X, A*) = R(X, A)* 
for all A in the resolvent set of A (and A*). Therefore, since A is a sectorial operator, 
being the generator of a bounded analytic semigroup, it follows that A* is a sectorial 
operator as well, and hence A* generates a bounded analytic semigroup (not neces- 
sarily strongly continuous at 0) on C^fOjOo) = Mf,[0,oo) by [2j Theorem 3.7.1] which 
coincides with T t * for all t > by the uniqueness of the Laplace transform. There- 
fore, the restriction of A* to D(A*) = L^C^oo) (i.e., the operator A°) generates a 
strongly continuous bounded analytic semigroup T® on L 1 [0, oo) by [2J Remark 3.7.13]. 
Furthermore, [H Remark 3.7.20] shows that for all n G IN and t > we have that 
Till G D((A*) n ) for all /i G M b [0,oo) and that If / G D((A ) n ) for all / G L^oo). 
Since D(A*) C L^O, oo), it follows that T> G L^O, oo) for all s > and /i G M 6 [0, oo). 
Therefore for allt > and \i G M&[0, oo) we have 

17 H = TtTtn = TfTtn e D((A ) n ) 

2 2 2 2 

for all n G IN. Thus by taking n = 8 x \t follows that T*^ G D((A°) n ) for all n. Using 
Theorem WM we have that D(A & ) = {g G L^oo) : D« y g(y) G ^[0, oo), D°~^(0) = 
0}, and then it follows that )" T t *5 x G L^C^oo) for all n. In particular, the 
transition probability distribution T*5 X has a density function y H- y, t) for all 
t > and all x ^ 0. To see that this function is smooth note that any positive integer 
m can be written in the form m = not — (3 for some integer n ^ 1 and some positive real 
number /3 < a. A straightforward calculation shows that for / G D((A & ) n ) we have 
[D°L y ) n f = D™f and I^ y D n _ a y f = D™" 13 f G L^O, oo) for all /? < a and n > 1. Since 
£>™ / = (-l) m (c//c/y) m /, we have {d/dy) m p{x, y, t) G #[0, oo) for all m. □ 
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Corollary 4.8. Let {/i n } C C®[0, oo) = M ac [0, oo) such that /i n — > n vaguely as 
n — > oo for some ji G M&[0, oo), then T®/i n — > T^\i vaguely as n — > oo. 

Proof. Since T®/x n = T t */x n for all n ^ 1, for all G Coo[0, oo) we have 



<j)(x)[Tf n n ](dx) = / 0(x)[T i */i n ](rfx) = [T t (f)](x)[i n (x)dx 



n— ^oo 



> \T t <l>]{x)n{dx) = U{x)[T^}{dx) 



(4.29) 



and the result follows. □ 

In Section we will apply Corollary 14.81 to compute the transition densities of the 
reflected stable process to any desired degree of accuracy, by solving the forward equa- 
tion numerically. For any initial state Z s = x, we will approximate the initial condition 
fio = 5 X in the numerical method by a sequence of measures \i n with L 1 -densities, and 
then Corollary 14.81 guarantees that the resulting solutions converge to the transition 
density in the supremum norm as n — > oo. 

Remark 4.9. The calculation (I4.26P of the forward generator (14.231) uses the L 2 [0,oo) 
adjoint formula 



f{x)A*g{x) dx = Af(x)g(x) dx. 
'o Jo 
The negative Riemann-Liouville fractional derivative A* = D°_ x is the L 2 (R) adjoint 
of the positive Riemann-Liouville fractional derivative A = D% (e.g., see [37]), but the 
adjoint calculation has to be modified on the half-line, using the boundary condition at 
the origin. This boundary condition is what distinguishes the forward equation (I4.28P 
of the reflected stable process Z t in (13.31) from the forward equation (14.71) of the original 
stable process Y t , since both have the same forward generator A* = D°_ x . This is a 
natural extension of the special case a = 2, where Z t is a reflected Brownian motion. 

Remark 4.10. Using integration by parts, one can write the backward generator in the 
form of an integro-differential operator (e.g., see Jacob 



(4.30) Af(x) = b(x)f'(x) + I [f(x + y)- f{x) - yf{x)\ <j>(x, dy), 

with coefficients 



x i-a 



(4.31) 

(/>(x,dy) = r(2_ a ) M 1 "dyh-Bfliv) + r(2 - a) X a£ - x ^ dy ^ 

The jump intensity (f>(x,dy) describes the behavior of the process Z t , which truncates 
jumps of the stable process Y t starting at the point x > in the state space, so that a 
jump (they are all negative) of size \y\ > x is changed to a jump of size x. This keeps 
the sample paths of Z t inside the half-line [0, oo). Since the drift b(x) is unbounded, the 
existence of a Markov process Z t with generator (I4.30p would not follow from general 
theory (e.g., see [THl Section 4.5] or [321 Section 3]). Hence the reflected stable process 
Z t is an interesting example of a Markov process with unbounded drift coefficient. 
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5. Transition density of the reflected stable process 



For the reflected stable process Z t started at Zq = 0, the form of the transition 
density was established in Proposition 14.11 To the best of our knowledge, there is 
no known analytical formula for the transition density y H- p(x; t, y) of the reflected 
stable process y = Z t+S started at x = Z s > 0. In this section, we compute and 
plot this transition density, by numerically solving the associated forward equation 
(I4.28p . The existence and smoothness of these transition densities is guaranteed by 
Corollary 14.71 Corollary 14.81 shows that T®g n (y)dy — > T^S x (dy) in the supremum norm 
as n — )• oo for any sequence of functions g n G L^C^oo) such that g n (y)dy —> 5 x (y)dy 
(vague convergence). We take g n {y) = n l[ x , x +i/n](y)- Then the solutions T®g n provide 
estimates of the transition density to any desired degree of accuracy. Theorem 14.51 shows 
that T®g n solves the forward equation (14.281) . Hence, we can compute the transition 
densities of the stable process by solving the forward equation numerically, with this 
initial condition. 

In order to compute the probability density p(x, y, t) numerically, we consider the 
fractional boundary value problem 

(5.1) d t u(y,t) = D a _ y u(y,t); u(y,0) = 5 x (y); LTfu&t) = 0. 

We develop forward-stepping numerical solutions Uh{yi,t) that estimate u(yi,t) at loca- 
tions yi = ih for i = 0, 1, . . . , N over an interval [0, y max \ in the state space, where y ma x is 
chosen large enough and h is chosen small enough so that enlarging the domain further, 
or making the step size smaller, has no appreciable effect on the computed solutions 
(e.g., the resulting graph does not visibly change). We approximate the delta function 
initial condition u(y, 0) = 5 x (y) by setting Uh(yi,0) = 1/h for yi = x and Uh(yi, 0) = 
otherwise, a numerical representation of the initial condition g n (y) = n t[ x , x +i/ n ](y) 
with h — l/n. 

Numerical methods for fractional differential equations are an active area of research. 
One important finding [231 121] is that a shifted version of the Griinwald finite difference 
formula ( 15. 6 p for the fractional derivative is required to obtain a stable, convergent 
method. Hence we approximate 

^ N+l-i 

(5.2) D a _ y u{y u t) n — ^ w^u h (y i+k ^,t) where < := (-l) fc 

for i ^ 2. Note that this approximation of D° y u(yi,t) does not depend on the value of 
Uh at the boundary y = 0. We enforce the boundary condition at yo = at each step; 
i.e., we set 

JV 

Uh(yo,t) = -^2w%~ l u h (y k ,t) 

k=l 

so that 

N 

D*?u{y Q ,t) w 5>rV(?M) = 

fc=0 
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since w, 



a-l 

o 



1. Finally, for i = 1 we approximate 



1 N 



fc=0 



(5.3) 



1 / - * \ 

= 7^ ( z2 w k u h(yk,t) -J2 w k lu h(yk,t) \ 

\k=l k=l / 



1 N 

-^J2 w k-l u h(yk^) 



k=l 



a-l 



using an elementary identity for fractional binomial coefficients w a — 

This leads to the following linear system of ordinary differential equations, 

fu h (y l7 t)\ (-\ a -i (u h {yut)\ 



-w 



a-l 
k-1 ■ 



(5.4) 



d_ 
dt 



1 



(-1 a-l .. 

1 —a w\ 

1 



\u h (y N ,t)J 



w 



N-l 



W 



N-2 







-a 



7 



\u h (y N ,t)J 



whose solution can be approximated using any numerical ODE solver (indeed, this 
is a linear system of the form u' = Au, so it has a solution u{t) = e tA u(0) for any 
initial condition). The resulting numerical solutions with y max — 12, h — 0.01, starting 
points x = 0,1,2,4, and a = 1.2,1.8 are depicted in Figures [1] and (2], where each 
frame represents a snapshot at times t = 0.5, 1 and 2 respectively. The L 1 -error of 
the numerical solutions for x = 0, compared to numerical estimates based on (14.11) and 
standard numerical estimates of the stable density, decays linearly with h. For h = 0.01 
the L 1 -error is less than 0.05 for a = 1.2, and less than 0.004 for a = 1.8, for every 
case plotted. A short Matlab code to compute the numerical solution is included in the 
appendix. 



Remark 5.1. It is interesting to note that the matrix in (15.41) is essentially the rate 
matrix of a discrete state Markov process in continuous time. Extending the state 
space to iV = oo, we obtain a Markov process on the state space {ih : i > 0} that 
approximates the reflected stable process, with Uh(xi,t) = P(Zf = ih). The transition 
rate from state ih to state jh for i > j > 1 is wf_j +1 ~ a(a — — j) _a_1 /T(2 — a), 
the jump intensity of the stable process Y t , in view of [26l Eq. (2.5)]. The transition 
rate from state ih for % > 1 to state jh for j = 1, in the first row of the rate matrix, is 
wfZi ~ i~ a /T(l —a), the rate at which the process Y t would jump into the negative half- 
line. This can be computed as <p{— oo, —ih) where <f)(dy) = a(a — l)\y\~ l ~ a dy /T{2 — a) 
is the Levy measure of the process Y t , e.g., see [261 Proposition 3.12]. 



Remark 5.2. The fractional boundary condition in (15. ip is a natural extension of the 
boundary condition for Brownian motion on the half-line. Let u(y,t) denote the 
transition density of reflected Brownian motion, the process Z t when a = 2. Then 
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Figure 1. The transition densities of p(x, y, t) with a = 1.2, x = 0, 1, 2, 4 
(left to right) at times t = 0.5 (left panel), t = 1 (middle panel), and t = 2 
(right panel). 




Figure 2. The transition densities of p(x, y, t) with a = 1.8, x = 0, 1, 2, 4 
(left to right) at times t = 0.5 (left panel), t = 1 (middle panel), and t = 2 
(right panel). 



d t u(y,t) = dyU(y,t) and 
(5.5) 



u(y + hA) — u(y,t) 
lim V ^ - ; = at x = 0, for all t > 0, 
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see for example Ito and McKean [THl Eq. 8]. The fractional boundary condition in ( 15. II) 
can be written in Griinwald finite difference form [261 Proposition 2.1] 



using the (fractional) binomial coefficients defined in (15. 2ft . When a = 2 we have 
= 1, u>" _1 = — 1, and u^ -1 = for A; > 1, so that ( 15.61) reduces to the classical 
condition ( 15. 5p . i.e., the one-sided first derivative. In either case (a = 2 or 1 < a < 2), 
the boundary term enforces a no-flux condition at the point y = in the state space. 

Remark 5.3. Bernyk, Dalang and Peskir [9j Appendix] computed the backward gener- 
ator of a general reflected stable Levy process. Caballero and Chaumont [TTJ Theorem 
3] compute the backward generator of a killed stable Levy process. It may be possible 
to develop the forward equation and compute the transition density for those process, 
using the methods of this paper. This would be interesting for applications to fractional 
diffusion, since it could elucidate the relevant fractional boundary conditions. 



The authors would like to thank Zhen-Qing Chen, University of Washington, Pierre 
Patie, Universite Libre de Bruxelles, and Victor Rivero, Centro de Investigation en 
Matematicas, for helpful discussions. 



The following Matlab code computes the transition density p(x, y, t) in the case x = 2 
for the reflected stable process Z t defined by (13. 3p . where Y t is a stable Levy process 
with characteristic function (13.21) and index 1 < a < 2. This code was used to generate 
the plots in Figures [T] and [2J 

/o°/„°/> Matlab script to compute p(x,y,t) 
°/ °/„ enter variables 

alpha=1.2; ymax=12; N=1200; t=[0, .5,1,2] ; x=2; 
°/o% initialise parameters 

h=ymax/N ; y= (h : h : ymax) ' ; 

uO=zeros(N, 1) ;uO(f loor(x/h)+l)=l/h; % initial condition 
°/ °/ Make Grunwald matrix 
w=ones(l ,N+1) ; 
for k=l:N 

w (k+1 ) =w (k) * (k-alpha- 1 ) /k ; 
end 

w=w/h~ alpha; 

M=spdiags (repmat (w,N, 1) , -1 : 1 : N-l ,N,N) ; °/ enter w's along diagonals 
M(l, : )=-cumsum(w(l :N) ) ' ; %change first row for BC 
%% Solve ODE system 

[~,p]=odell3(@(t,u) M*u,t,uO); 



(5.6) 




at y — 0, for all t > 0, 



fc=0 
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